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Simulations of minor mergers. I. General properties of 
thick disks 
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ABSTRACT 

We present simulations of the formation of thick disks via the accretion of two- 
component satellites onto a pre-existing thin disk. Our goal is to establish the detailed 
characteristics of the thick disks obtained in this way, as well as their dependence on 
the initial orbital and internal properties of the accreted objects. We find that merg- 
ers with 10-20% mass of the mass of the host lead to the formation of thick disks 
whose characteristics are similar, both in morphology as in kinematics, to those ob- 
served. Despite the relatively large mass ratios, the host disks are not fully destroyed 
by the infalling satellites: a remaining kincmatically cold and thin component con- 
taining ~ 15 — 25% of the mass can be identified, which is embedded in a hotter 
and thicker disk. This may for example, explain the existence of a very old thin disk 
stars in the Milky Way. The final scale-heights of the disks depend both on the initial 
inclination and properties of the merger, but the fraction of satellite stellar particles 
at ~ 4 scale-heights directly measures the mass ratio between the satellite and host 
galaxy. Our thick disks typically show boxy isophotes at very low surface brightness 
levels (> 6 magnitudes below their peak value). Kinematically, the velocity ellipsoids 
of the simulated thick disks are similar to that of the Galactic thick disk at the solar 
radius. The trend of ozI<Jr with radius is found to be a very good discriminant of the 
initial inclination of the accreted satellite. In the Milky Way, the possible existence of 
a vertical gradient in the rotational velocity of the thick disk as well as the observed 
value of az /<Jr a t the solar vicinity appear to favour the formation of the thick disk 
by a merger with either low or intermediate orbital inclination. 

Key words: methods: numerical - galaxies: formation, kinematics and dynamics, 
structure - Galaxy: disc, kinematics and dynamics, formation 



1 INTRODUCTION 

The different components of disk galaxies contain key in- 
formation about various stages of the formation history of 
these systems. In this sense, one of the most significant com- 
ponents for studying signatures of galaxy formation are the 
thick disks because they contain imprints of the state of 
early disks and their interaction with th e galactic environ- 
ment l|Freeman fc Bland- Hawthorn! [200^1 . 

Until now thick-disks have been detected in SO 
galaxi es (pBursteinl [l979l: iTsikoudl 1 19791 ). in the Milky 
Way dGilmore fc Reidl 1 19851) and in many other s p rals 



llvan der Kruit fc Searld Il981allb1; Ijensen fc Thuanl Il982l: 
van Dokkum et al.lll994;lMorrison et al.ll 19971: iPohlen et al. 
2000; iDalcanton fc Bernsteinl 120021 : lYoachim fc Dalcanton 
20061 . and references therein). As such, thick disks appear 



to be a rather ubiquitous structural component of galaxies. 
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Historically, there have been two types of scenarii pro- 
posed to explain the formation of thick disks: (i) "dissipa- 
tional", and (ii) "predominantly dissipationless" models. In 
the first case, thick-disk stars are formed during the dissi- 
pational collapse of gas with a large scale-height, after the 
halo has formed and before the thin disk has completely de- 
veloped. A variant of this model, more in line with modern 
cosmolog y, is given bv lBrook et"al] i|2004l2005l ) and consists 
of the formation of the thick disk in an epoch of multiple 
mergers of gas-rich building blocks. In this class of models, 
one may expect a smoo th transition between properties of 
the thin and thick disks l|Eggen et~a l, 1962; G ilmore fc Wvse 



1986; 



1995; 



Norris fc Ryan 199ll: iBurkert et al.|[l992l : iPardi et~ 



Fuhrmannl 2004 ) . Evidence of such a process at work 



may be th e chain and clumpy galaxies o bserved at high- 
redshift bv lElmegreen fc Elmegreenl l|2006h . which have also 
been suggested to be the progenitors of thick disks (e.g. 
iBournaud et al.ll2007al ). In the "predominantly dissipation- 
less" models the thick disk stars were either: (1) vertically 
"heated" from a pre-existent thin disk during a (signifi- 
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cant) minor merge r llQuinn et ah 1 1 19931; iMihos et al J 1 19951; 
Walker et all Il99rj; [Robin et all 1 19961; I Velazquez fc White 



199S : lAguerri et al.ll200ll ; IChen fc the SPSS Collaboration 
200 if ); (2) directly deposited at large scale-heights as 



tidally stripped debris during the accretion of smaller 



satellite galaxies (Bck ki fc Chiba|l200ll ; lGilmore et ajj 
Abadi et all 120031; iMartin et al.l |2004|; iNavarro et al 



2002 



2004 



Hclmi et al.l 2006); or (3) the product of the dissolution of 



massive thin-dis k clusters wit h small radii and large veloc- 
ity dispersions (|Kroupal l2002). The first and third models 
have in common the presence of a pre-existing disk compo- 
nent. In this second class of models, the thick disk may be 
characterised as a completely foreign component. 

Models involving either one or more minor mergers - 
including the case of gas-rich accretion proposed by Brook 
et al.- are natural in the context o f current theories of hier - 
archical structure formation (e.g., Kazantzidis et al"1l2007l ). 
Such models are also supported by modern studies of both 
kinematical and chemical properties of the thick disks in the 
Milky Way and in external galaxies (lYoachim fc Dalcantonl 
120051 . 120061 ; ISeth et al.ll2005l ; lMouldll2005l ). The observed lack 
of vertical colour and chemical gradients, as well as the pres- 
ence of counter-rotating disks, favours scenarios in which 
thin and thick disks formed as separate entities. However, 
it is unclear to what extent such models reproduce the de- 
tailed properties of observed thick disks. For example, in the 
dissipationless minor merger scenario, the thick disk should 
contain stars from both the heated thin disk and the ac- 
creted satellite. What fraction of the stars come from each 
constituent? If they are mostly satellite debris, one may ex- 
pect a significantly metal-poor thick disk, which appears 
to be i ncons istent with the observations reported in, e.g., 
iMouldl (|2005l ). 

It is the lack of very detailed predictions that motivates 
us to revisit the problem of thick disk formation. In this 
paper, we explore the hypothesis of the heating of a pre- 
existing thin disk by a single minor merger. In the scenario 
we envisage a new thin disk component would form from 
the accretion of cold gas after the merger has taken place. 
However, here we only study the global properties of the fi- 
nal merger product. In a follow-up paper we will focus on 
its phase-space structure with the aim of making a detailed 
comparison to the thick disk in the Solar neighbourhood, 
and eventually to uncover the debris from the object that 
may have triggered its formation (Villalobos & Helmi, in 
prep). It is important to note that we shall neglect the ef- 
fect of gas on the properties of the remnant thick disk for 
the moment. This simplification should be borne in mind, 
in particular since the structure of minor merg e r rem nants 
may be different, as suggested by iNaab et al. (120061 ) and 
Younger et all (|2008l ). 

Significant work has been carried out in the past to 
model the disk heating process by a mi nor merger from the 
numer ical point of view, starting from Quinn fc Goodman! 
1986) (QG86I); iQuinn. Hernquist. fc Fullagarl (Il993fl 



QHF93h : 



Mihos et al 
— I — 



(ll995l);IWalker. Mihos. fc Hernquistl 
19961) (IWMH96I); iHuang fc Carlberd (Il997l^ (IHC97I) ; 
ISellwood et all dl998h7fVerazQuez fc White! (Il999l ) l|VW99l ). 
and more recently iKazantzidis et al.l (I2007T 1. and simulta- 
neously with this work iRead et alj (|2008j ). These papers 
have essentially shown that it is relatively easy to produce 
thick disks whose general properties are consistent with 



those observed. However, we believe there is still room for 
improvement, both in the initial conditions used to model 
this process, as well as in the degree of detail necessary for 
comparisons to the latest observations of thick disks. 

Besides the work mentioned above, in recent years sev- 
eral authors have studied minor mergers with a number of 
different goals (other than the formation of thick disks). 
These include for example, studies of the merger remnants 
produced by encounters between disk gala xies (as a way of 
prod u cing a bulge dominated system, e.g 



120031 ; IJesseit et al] 120051 ; INaab fc TruUuol 



Naab fc Burkertl 
2006). More 



cently, also the effects of gas, star formation and feed- 
back have been t aken into account dBournaud et~al Il2005l ; 
INaab et al l 120061 ; Idi Matteo et~ail 120071 ; ICox et all (20081 ) . 
The eff ects of repeated minor mergers have been consider ed 
in, e.g. iBournaud et al. I (|2007bl) ; IKazantzidis et ail (|2007l ). 

The space of initial conditions to tackle the formation 
of the thick disk via minor mergers is very large, and it is not 
desirable to probe it randomly. We have therefore made the 
following physically motivated choices, (i) We model the for- 
mation of thick disks at two different redshifts, by scaling the 
properties of the host galaxy an d accreted satellite accord - 
ing to cosmological models, as in lMo. Mao, fc White! (|l998l ) 
(ii) We consider the accretion of relatively massive satellites 
(10% or 20% mass ratios), embedded in dark-matter halos, 
and with stellar distributions that are initially either spher- 
ical (and on the fundamental plane of dE+dSphs galaxies) 
or disky. (iii) The satellites are released much farther away 
from the host disk compared to previous studies, and their 
orbits are consistent with thos e of infalling s ubstructures in 
cosmological simulations (e.g., Benson 2005). 

The outline of this paper is the following: In §2 we 
describe in detail the numerical procedure used to build 
self-consistently the different components of both the host 
galaxy and satellites including our choices for the numerical 
parameters and the orbital parameters of the satellites. §3 
describes the outcome of the simulations, focusing on the fi- 
nal properties of thick disks and the evolution of satellites. In 
§4 we present a discussion and limitations of our approach. 
In 85 we summarise the main results. 



2 SETTING UP THE SIMULATIONS 

In this section, we describe in broad terms the proce- 
dure adopted to model the host disk galaxy and the (to 
be accreted) satellite. We refer the interested reader to 
the Appendix for more details. We consider two configu- 
rations: a merger with a host whose properties resemble the 
Milky Way today (our "z=0" experiment), and a merger 
with a smaller host disk galaxy (this is our "z—1" experi- 
ment). The "z=0" configu ration has been often used in the 
past in the same cont e xt dQG86l: IQHF93I; IWMH96I; l"HC97l; 
VW99I; lAguerri et al.l I2OO1I ; iFont et all l200ll ; lArdi et all 
20031 ; lHavashi fc Chibal 120061 ) . In the "z=l" configuration 
both the host system and the satellite's properties are scaled 
to tho se expected at z=l according to the model of lMo et alj 
(1998). In this case, the aim is to simulate the merger event 
that might have given rise to the Milky Way thick disk. In 
this configuration the mass of the present-day thick disk of 
the Milky Way is roughly equal to the combined mass of the 
host disk galaxy and the stellar component of the satellite. 



Table 1. Properties of host disk galaxies. 
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Table 2. Properties of satellite galaxies. 
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"z=0" "z=l" 



NFW Halo 


Virial mass, M v i r 


10 12 


5.07 x 10 11 


[Mq] 


Virial radius, R v i r 


258.91 


122.22 


[kpc] 


Concentration, c 


13.12 


6.56 




Circular velocity, V c (R v i r ) 


129.17 


133.87 


[km/s] 


N H 


500000 


500000 




Softening, t halo 


0.35 


0.41 


[kpc] 


Disk 


Disk mass, M<ji s k 


2.8 x 10 10 


1.2 x 10 10 


[Mq] 


Scale-length, Ftp 


3.5 


1.65 


[kpc] 


Scale-height, zq 


0.35 


0.165 


[kpc] 


Toomre Q(R= 2AR D ) 


2.0 


2.0 




N D 


100000 


100000 




Softening, t disk 


0.05 


0.012 


[kpc] 



In this section we also describe the orbital parameters 
of the various experiments. Furthermore, we explain the 
choices made for the numerical parameters employed in our 
simulations, and describe the global stability of the system. 

2.1 Main disk galaxy 

We model the main disk galaxy as a self-consistent two- 
component system, containing a NFW dark matter halo 
(N avarro et al .1119971 ) and a stellar disk. The dark halo is adi- 
abatically contracted in r esponse to the formation of a stel- 
lar di sk in its central part l|Blumenthal et aTl il 1986; M o et all 
1 19981 ). The disk co mponen t is constructed follow ing the pro- 
cedure outlined by QH F^9ll and lHernquistl (| 19931 ). and follows 
a density profile of the form: 

pd{R ' z) = sJp^v^tI) sech2 (i) (1) 

where Md is the disk mass, Rd is the exponential scale- 
length, and zp is the exponent ial scale-height. 

Following IMq et all (|l998f l. the ratios M haio /M dlak and 
Rvxt/Rd have been kept nearly constant for all redshifts. 
The scale-height zo has been kep t as O.lRp accord ing to 
observations of external galaxies l|Kregel et al.l l2002). This 
has also been assumed in our "z—1" experiments. By keeping 
this ratios constant, the way in which the mass and dimen- 
sions of the disk component scale with redshift is simplified, 
since it follows the same scaling with redshift as the halo 
within which it is embedded. The chosen values for the pa- 
rameters of the main disk galaxy are listed in Table [T] for 
our "z=0" and "z=l" experiments. 

2.2 Satellite Galaxies 

The satellite galaxies are designed self-consistently with 
both dark matter and stellar components. The dark halo 
follows a NFW density profile, and the initialisation proce- 
dure is identical to that for the host system. 

1 Note that sech 2 (z/2zo) ~ ex p( — l z l/ z o) f° r l 2 l 3> zo- 





"z=0" 


"z=l" 




NFW Halo 


Virial mass, M v i r 


2 X 10 11 


1.01 x 10 11 


[Mq] 


Virial radius, R v i r 


151.40 


71.35 


[kpc] 


Concentration, c 


16.18 


8.09 




Circular velocity, V c (R v ir) 


75.50 


78.10 


[km / s] 


N H 


100000 


100000 




Softening, t halo 


0.14 


0.07/0. 14 1 


[kpc] 


Disk 


Disk mass, M^isk 


5.6 x 10 9 


2.4 x 10 9 


[Mq] 


Scale-length, Rd 


1.69 


0.96 


[kpc] 


Scale-height, zq 


0.17 


0.095 


[kpc] 


Toomre Q(R = 2AR D ) 


2.0 


2.0 




N D 


100000 


100000 




Softening, e disk 


0.024 


0.007 


[kpc] 


Bulge 


Bulge mass, M bu i ge 


5.6 x 10 9 


2.4 x 10 9 


[Mq] 


Scale radius, at, 


0.9 


0.709 


[kpc] 


Velocity dispersion, (Jq 


96.29 


69.06 


[km/s] 


N B 


100000 


100000 




Softening, e bulge 


0.07 


0.07 


[kpc] 



; Softenings used in disky and spherical satellites respectively. 



We consider two possible stellar distributions: an expo- 
nential disk and a sph erical Hernquist bulge. In this case, 
the density is given by (|Hernquistlll990h : 



where Mb is the bulge mass and at is the scale radius. 

Both types of satellites satisfy M to tai,sat = 
0.2M tota i,host for the "z=0" and "z=l" experiments. 
In particular, the mass ratio between dark matter and 
luminous matter in the satellites is set to be the same as in 
the host galaxy. This implies that the mass of the stellar 
components of our satellites is 20% of Mdisk,host, similar to 
those ado pted in previous w o rk on disk hea ting by satellite 
accretion l|QHF93l; |WMH96| ; IHC97I : IVW99T ) . Note however 
that the initial total mass of the satellite is comparable to 
that of the host disk. For completeness, we also include the 
case of a satellite with a stellar mass of 10% of Mdisk, host- 
In this case the total mass of the satellite is also 10% of 
that of the host galaxy. 

The spherical satellites li e on the observed fun damental 
plane of dE+dSphs galaxies l|de Riicke et al.ll2005t ): 

log L B ~ 4.39 + 2.55 log a (3) 

log L B ~ 8.65 + 3.55 log R e (4) 

where L B is the blue-band luminosity, and <to is the cen- 
tral velocity dispersion. The effective radius is related to 
the scale radius by R e « 1.82at, for a Hernquist density pro- 
file. Finally, we fix the mass-to-light ratio T B = 2Ts to 
derive the stellar mass of our spherical satellites. 

The structural parameters of the disky satellites, Rd 
and zo are set in the same way as for the host. 

The above relations fully specify the properties of our 
satellites, both for the "z=0" and "z=l" experiments. Note 
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that the variation with redshift is always linked to that of 
the host halo. Table [2] lists the properties of both spherical 
and disky satellites for "z=0" and "z=l". 



2.3 Orbital Parameters 

We release our satellites f rom significantly larger distances 
than previous works (e.g.. lQG86l ; IVW99l ; iKazantzidis et al.l 
120071 ) . For the "z=0" case, the satellite is launched from 
the virial radius of the host galaxy computed at z=l 
(R vir = 122.22 kpc » 35R D )- For the "z=l" experiment, 
it is launched from a distance of 83.9 kpc (w 50Rd) which 
corresponds to the virial radius of the host galaxy computed 
at z—1.6. 

T o initiali s e the orbital velocities of the satellites, we 
follow iBensonl (|2005l) . Benson determined the orbital pa- 
rameters of DM substructures at the time they crossed 
the virial radius of their h ost halos (see also iTormenl 1 1997*1 ; 
iKhochfar fc Burkertl [2006" ) ■ We choose for the velocities of 
our satellites the most probable values of the distributions 
as given for z=0 (since little variation is visible as a function 
of redshift). The radial and tangential velocity distributions 
peak, respectively, at 0.9 and 0.6 in units of V c (R v i r ). 

We consider for our satellites initial orbital inclinations 
of 0°, 30° and 60° (with respect to the plane of the host 
disk), in both prograde and retrograde directions with re- 
spect the rotation of the host disk. In the case of disky 
satellites, their midplanes are parallel to that of the pri- 
mary disk, implying that they are, respectively, inclined by 
0°, 60° and 30° with respect to their orbital planes. Table [3] 
summarises the orbital parameters adopted for both spheri- 
cal and disky satellites, for the configurations at "z=0" , and 
"z=l". 

Note that within each configuration, all the satellites 
are initially released from the same distance with the same 
velocity modulus, i.e., they all have the same total energy 
and modulus of the angular momentum, but the latter vary- 
ing its orientation. This implies that the initial apocentre, 
pericentre (and hence orbital eccentricity) are the same for 
all experiments. However, as we shall see in £13. li the orbits 
evolve due to dynamical friction, with some dependence on 
the internal properties of the satellite. Therefore, at the time 
the satellite merges with the disk, the orbits of all our ex- 
periments are different. 

With these initial conditions, we find that our satellites 
have completely merged with the host disks by z~0.4 and 
by zRil, respectively for the "z=0" and "z=l" cases. 



Table 3. Initial orbital parameters of satellite galaxies. 



2.4 Numerical Parameters 



The i V-body systems are evolved using Gadget-2.0 (|Springell 
120051 ) a well documented massively parallel TreeSPH code. 
Our choices of the numerical parameters (number of par- 
ticles iV in each component in the system; softening e and 
timestep At) are described in detail in £|A2I Tables [T] and [2] 
list the values used for each component in our simulations. 
The maximum timestep (not listed in the Tables) is 0.25 
Myr. Typically the energy and angular momentum are con- 
served to better than 0.1% over 9 Gyr of evolution for our 
main disk galaxy configured at "z=0" . 



i 


X 


z 


v x 






"z=0" 


0° 


122.2 


0.0 


-137.7 


0.0 


11219.8 


.30° 


105.8 


61.1 


-119.2 


-68.8 


9716.6 


60° 


61.1 


105.8 


-68.8 


-119.2 


5609.9 


"z=l" 


0° 


83.9 


0.0 


-118.2 


0.0 


6615.9 


30° 


72.7 


42.0 


-102.4 


-59.1 


5729.5 


60° 


42.0 


72.7 


-59.1 


-102.4 


3307.9 



NOTES: 

- Distances in kpc, velocities in km/s, angular momentum in 
kpcxkm/s. 

- In all cases, y = kpc, and v y = 91.8 km/s for "z=0" and 
v y = 78.8 km/s for "z=l". 

- Listed v y and L z are for prograde orbits. Retrograde orbits 
have the opposite sign. 

- Initially r apo = 77 kpc and 49 kpc, r per i ~ 10 kpc and 5 kpc, 
respectively for the "z=0" and "z=l" configurations (as mea- 
sured from the first apocentre and the subsequent pericentre). 
The corresponding eccentricities e = (r apo — r per i)/ (r apo + r per i) 
are 0.77 and 0.82. 



2.5 Evolution of Isolated Host Galaxy 

Before including the satellite, the host galaxy is simulated 
in isolation to test its stability in the absence of any external 
perturbation. As we show in ijA3lour host galaxies are ver; 
stable in their properties (cf., VW99I ; iGauthier et al.ll20pj " 
for the amount of time needed to complete the experiments. 
Therefore, we are now ready to focus on how these systems 
evolve when they suffer a minor merger. 



3 RESULTS 

In total, 25 simulations have been carried out to study the 
formation and global properties of thick disks as a result 
of the merger between a host disk galaxy and a satellite. 
The simulations explore combinations of the following ele- 
ments: two configurations for the progenitors ("z=0" and 
"z=l"); two morphologies for the stellar component of the 
satellite (spherical and disky) ; two total mass ratios between 
the satellite and the host galaxy (10% and 20%); and three 
initial orbital inclinations for the satellite with respect to 
the midplane of the host disk (0°, 30° and 60°), in both 
prograde and retrograde directions. 



3.1 Orbital Evolution of the Satellites 

To study the evolution of the orbits in our experiments we 
follow the location of the centre of mass of the satellite (CM) 
with respect to the host galaxy. The CM is defined as the 
mean position of bound stellar and DM particles of the satel- 
lite. Fig. [T] shows the trajectories of the CMs for the ex- 
periments configured at "z=0" and "z=l" with a spherical 
satellite on a prograde orbit with initial inclination i = 30° . 
The XY projections clearly illustrate how radial the orbits 
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Figure 1. XY and XZ projections of the trajectory of the centre 
of mass of a spherical satellite as it decays towards the centre of 
the host galaxy in our "z=0" and "z=l" experiments with initial 
inclination 30°. The coordinate system is centred on the centre 
of mass of the host disk. The inset panels show in more detail the 
trajectory of the satellite at late times before it is fully disrupted. 
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Figure 2. Evolution of the radial separation between the centres 
of mass of the host disk and the satellite, and of the z-distance 
of the satellite's centre of mass with respect to the host disk 
plane, for spherical satellites on prograde (top panels) and retro- 
grade (bottom panels) encounters. The various curves correspond 
to different orbital inclinations for the "z=l" configuration, and 
the 10% experiment is shown for the inclination of 30°. The in- 
sets show in more detail the trajectory of the satellites at late 
times. The lack of a clear trend in the amplitude of the verti- 
cal oscillations is likely due to the difficulty in determining the 
exact orientation of the disk plane at the ~ 1 kpc level. Similar 
behaviours are observable for the "z=0" case. 



are and how rapidly they decrease in amplitude due to dy- 
namical friction between the satellite and the host system. 
Note that these orbits are quite different from the circu- 
lar ones u sually used in earlier studies (e.g., QG86; QHF93; 
IWMH96I ). The XY and XZ projections show that the satel- 
lites stay on their original orbital planes as they decay, until 
they enter into the zone dominated by the host disk, when 
they are drawn onto the disk plane and spiral in towards 
its centre. This is particularly clear for the lighter satellite 
in our "z=l" experiment that spends ~1 Gyr on the disk 
plane before sinking in further. 

Note that most of the angular momentum of the "satel- 
lite + host system" is in the orbital motion of the satellites. 
This is because the satellites are relatively massive and have 
a very large initial distance from the host. As a result, when 
the satellite decays in an inclined orbit the disk is strongly 
tilted in both the prograde and retrograde cases (see §3.3). 

As expected, the trajectory of the lighter satellite is 
more extend ed because dynamical fri ction is less efficient 
in this case ijBinnev fc Tremainei [19871 ). Our 20% satellites 
decay completely after ~3 Gyr in the case of "z=0" experi- 
ments, and after ~2 Gyr for "z=l" experiments. In compar- 
ison, the satellite with half of the mass takes the double of 
time to sink starting from the same initial orbital parame- 
ters. 

Initially, the orbital decay is due to the dynamical fric- 
tion against the host halo. This implies that the decay rate 
does not depend on inclination, orbital direction or stellar 
mass distribution (see Fig. [2]). On the other hand, when the 
satellites approach the centre of the system (where the disk 
is a significant contributor to the global force field) , the orbit 
decays by dynamical friction also against the host disk. Thus 
at later times, prograde low inclination or bits decay faster 
than re trograd e or high inclination orbits (jOG86t IWMH9fj| ; 
see also lHC97T ). 

3.2 Satellite Mass Loss 

Fig. [3] shows the mass loss evolution of both DM and stellar 
components of spherical and disky satellites in the "z=0" 
and "z=l" experiments. 

In order to calculate how much mass remains bound to 
the satellite at a given t ime we implemente d in Gadget-2.0 
the following procedure (|Benson et alj|2004 ): 

(i) Start by considering all the satellite particles that were 
bound to the satellite at the previous timestep (or simply 
all satellite particles for the first timestep). 

(ii) Compute the mass of the satellite from these particles 
along with the position and velocity of the CM. 

(iii) For each particle in this set, determine whether it is 
gravitationally bound to the other particles in the set. 

(iv) Retain only those particles that are bound and go 
back to step (ii). Repeat until the mass of the satellite has 
converged. 

Since the satellites were initialised in the absence of an 
external potential, as soon as they are placed within the 
host potential a large fraction (70%) of the more extended 
DM component rapidly becomes unbound before the first 
pericentric passage. After that, the mass loss rate of the 
DM component mostly depends on its initial mass. 

The mass loss rate of the stellar components depends 



6 A. Villalobos and A. Helmi 



strongly on in itial mass, and orbital parameters l|WMH96l; 
HC97; VW99) but also on the stellar mass distribution. As 
satellites decay, prograde orbits with lower inclinations lose 
mass faster than retrograde orbits with higher inclinations, 
due to the stronger tidal interaction with the host galaxy. 
This trend is more notorious for disky satellites. Spherical 
satellites are also characterised by a more extended "knee" 
in the mass loss in comparison to disky satellites. The lighter 
satellite experiences a slower mass loss compared to heavier 
satellites, because it suffers less dynamical friction and hence 
is on a less bound orbit. 

Once a satellite has sunk onto the plane of the host disk, 
its fate will depend on its instantaneous mean density com- 
pared to the mean density of the host at a given location. If 
the mean density of the satellite is larger than that of the 
host system then the most bound particles of the satellite 
will reach the galactic centre as a distinctive core causing 
more damage to the host disk. Otherwise, the satellite will 
receive most of the damage, being heated and torn apart 
by the host disk. In our simulations only the mergers with 
spherical satellites in the "z=0" experiments deposit in the 
galactic centre final cores of up to 20% the initial stellar 
mass. These final cores are on the lighter side in co mparison 
with previous studies using high density satellites (|QHF93l . 
20%; IWMH9rj 45%). It is interesting to note that spherical 
satellites with higher inclinations give rise to the formation 
of less massive cores. This can be explained by the fact that 
satellites on higher inclinations experience more disk cross- 
ings through the host disk as they decay, compared to ones 
on lower inclinations. In this case disk shocks perturb the 
structure of the satellite and cause additional mass loss (see 
iBinnev fc Tremaine|[l987l ; l0HF93l ). 

Figure shows that in general most of the dark matter 
is stripped off early, and deposited at very large radii (see 
Fig. [2}. Therefore the fraction of dark matter accreted from 
the satellite and deposited in a disk-like struc ture is very 
small in our simulations, in comparison to what Rea d et all 
(2008) find. This may be explained by the fact that our 
satellites, although of comparable mass, are launched from 
much larger distances. 

3.3 Description of the Mergers 

Figure [4] illustrates the morphological changes in the host 
disk during merger events configured at "z=l". The initial 
inclination of the satellite is 30° for both prograde and ret- 
rograde orbits. 

As the satellite decays in its orbit, it induces the forma- 
tion of noticeable spiral arms in the host disk, which trans- 
port angular momentum from the central parts towards the 
outskirts. Once the satellite sinks onto the plane of the host 
disk it transfers kinetic energy from its orbit to the particles 
in the disk, increasing their vertical motions and causing a 
visible thickening. At the same time the disk responds to the 
decaying satellite, by tilting its plane in order to conserve 
the total angular momentum of the system (although a sig- 
nificant amount of the satellite's initial angular momentum 
has by this time, already been transferred to the host halo) . 

Figs. [5] and [6] show the distribution of stellar particles 
for both spherical and disky satellites on prograde orbits 
with initial inclination 30°. In the early stages of the orbital 
decay, the satellite is stripped leaving trails of particles on 
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Figure 3. Mass that remains bound for both the dark matter and 
stellar components of the satellites in our experiments at "z=0" 
(upper panels) and at "z=l" (bottom panels). The bound mass 
is in units of the initial mass of the corresponding component. 



orbits with inclinations similar to that of the satellite ini- 
tially. Since the stars are initially deeply embedded in the 
satellite's dark matter halo, only a small fraction of the stel- 
lar debris is deposited at large radii. Most of the stars from 
the satellite end up in a disk-like configuration, with the 
same orientation as the final disk, but one that is somewhat 
thicker and more extended (see §3.4). Noticeable shells of 
debris material are formed as time goes by. These struc- 
tures are a consequence of the i nteraction of a dynamicall y 
cold system with a larger one (|Hernquist fc Quinnl I1988T ). 
In general the survival of these shells will depend on the 
mean phase-space density of the infalling satellite and also 
on its orbit. In the case of spherical satellites shells are visible 
typically since t ~ 1.7 Gyr, lasting ~2 Gyr; and for disky 
satellites much sharper shells are seen starting at t ~ 1.5 
Gyr, being still noticeable by the end of the simulation, i.e., 
~2.5 Gyr later. Shells are rather common features related to 
merger events, being observed in many elliptical and spiral 
galaxies. An important characteristic of shells is that they 
usually survive for a long time in phys ical space, as previ- 
ous numeri cal studies have shown (e.g. iHernquist fc Quinnl 
1988, 1989). The presence of such a structures in the solar 
vicinity and the possible signatures imprinted on them dur- 
ing the formation of the thick disk will be explored in Paper 
II (note that such features have alr eady been propos ed to 
explain the Monoceros ring, see e.g. iHelmi et all [2003). 
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Figure 4. Evolution of the host disk (initially shown edge-on) during the merger with a spherical satellite with inclination of 30° for 
prograde and retrograde orbits in the "z=l" experiment. Only host disk particles are shown for clarity. Note that the disk appears over 
thickened or distorted because of projection effects, especially in the case of the retrograde orbit. The bottom panels show edge-on views 
of the final system. 



3.4 Properties of the merger products 

By the end of the simulations the morphological, structural 
and kinematical properties of the heated disks and satellite 
debris have settled down and do not evolve further. This 
occurs ~2 Gyr after either the satellite has been disrupted 
or its core has reached the centre of the host disk. This means 
that the properties of the final thick disks do not change after 
t=5 Gyr and t=4 Gyr for the systems configured at "z=0" 
and at "z=l" respectively, in the case of heavy satellites. 
This timescale is ~ 6 Gyr for the lighter satellite in our 
"z=l" experiment. 

We now study in more detail the characteristics of the 
final disks. To this end we use a reference frame centred on 
the centre of mass of the final product, and aligned with its 
principal axes in such a way that the rotation axis defines 
the ^-direction. 



3.4-1 Morphological properties 

The left panels in Fig.[7]show the morphologies of the heated 
host disks at the final time of the "z=l" experiments (t=4 
Gyr), for the case of the spherical satellites. Prograde orbits 



and lower inclinations induce on the host moderate arms and 
more radial expansion in comparison with retrograde orbits 
and higher inclinations. On the other hand, higher prograde 
inclinations are more efficient at thickening the disk, espe- 
cially in the outskirts. For instance, a prograde satellite with 
inclination of 60° only causes a slight increment of the radial 
extension compared to the coeval control disk, but induces a 
noticeable thickening compared to the same control disk and 
to the other inclinations. Satellites on retrograde orbits have 
a similar thickening effect on the disk but a considerably 
milder influence on the formation of tidal arms and radial 
expansion compared to satellites on low inclination prograde 
orbits. Notice also that the satellites durin g their decay can 
induc e the formation of weak bars (see also Berentzen et al.l 
l2004| y Some warping in the disks is also visible in the case 
of mergers with inclinatio ns of 60° for bo th prograde and 
retrograde orbits (see also lQHF93t IVW99I) . 

The panels on the right side of Fig. [7] show the thick 
disks obtained, now including the contribution of the satel- 
lite's stellar particles. Their final structure is dominated by 
the heated disk (compare to the left panels), except in the 
outer regions, where the contribution of satellite debris is 
important. The outskirts are clearly thicker for satellites on 



8 A. Villalobos and A. Helmi 



Prograde Face— on 



Prograde Edoe— on 



30 

ID 



-10 

-ao E- 

30 
10 


-10 

-ao F- 



M 1 11 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

0.50 Gyr 



"0.75 



r ■ | • ■ ■ ■ | t i^^flJ " " I 

1,25 ' 



30 

10 
D 

-10 
-20 



mm 




.1.. I. ...I.... I. 



1 1 1 11 1 1 



I 1 1 1 1 [ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

i.so: 




:f2.25 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
1.00 




■4.00 




iiiiiiui.iiii.il 



2 
10 


-10 
-30 

30 

10 


- 10 

30 



'I 1 11 1 I 1 1 1 1 1 1 1 1 1 I 1 1 1 1 I 1 

50 Gyr 



n |iiii|iiii|i N i|i n i|i )i i iiii|iiii|iiii| u ii|i :i |i h | n i |i h i|i h i|i ) H |ii N | n ii| nn | n ii| n] | iiii|i n i|iiii|iiii| ni i| i i h |iiii|iiii|iiii|i i 



30 
10 


-10 

-30 



|||I m ||II|||I m || 



1.35 



1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
i 1.50 




"2.00 




■0.75 



1 1 1 1 1 1 1 a ' 1 ■ 1 ■ ■ 1 ■ ■ 1 1 1 ■ 1 1 1 1 ■ 1 ■ 1 1 ■ 



--1.00 



II I ' | II I ' | I I I I | ' I I 

JT.75 



;f2.25 



1 ■ ... 1 ■ a: ■ 1 ■ ... 1 




:r4.oo 



-30 -10 10 20 -20 -L0 10 SO -30 -10 10 30 
x/kpc 



-30 -10 10 30 -30 -L0 10 30 -30 -10 10 30 
x/kpc 



Figure 5. Evolution of the spherical satellite with inclination 30° during the prograde merger with the host galaxy in the "z=l" 
experiment. In each snapshot the reference frame has been centred on the centre of mass of the host disk and has also been rotated to 
eliminate the tilting of the whole system with respect to the original frame. The labels face-on and edge-on are relative to the host disk 
(not shown here for clarity). 
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Figure 6. Same as Fig. [5] for the case of the disky satellite. 



higher inclination orbits, although their debris does not show 
the warp feature characteristic of the heated disks. Also no- 
ticeable is the difference in the distribution of satellite de- 
bris between prograde and retrograde orbits for the case of 
coplanar infall (in edge-on views). 



The contour levels shown in this figure have been drawn 
4.5, 6.2, 8 and 9.7 magnitudes below the central surface 
brightness of the remnant system. If we assume a mass- 
to-light ratio Ty = 2Tq,v for the host as well as for the 
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Figure 7. Face-on and edge-on views of the final morphologies of heated disks (left) and thick disks (right) at the end of the simulations 
in "z=l" experiments, 4 Gyr after the infall of the spherical satellite. In each case the tilting induced by the satellites has been eliminated 
to facilitate the comparison with the coeval control disk. The contours correspond to 4.5, 6.2, 8 and 9.7 magnitudes below the central 
surface brightness of the remnant system. For the experiments shown here, and assuming a mass-to-light ratio Yy = 2T© y , these would 
be located at 22.5 mag/arcsec 2 (innermost), 24.2 mag/arcsec 2 , 26 mag/arcsec 2 and 27.7 mag/arcsec 2 (outermost) in the V-band. 



satellite stars, these contours correspond to 22.5, 24.2, 26 
and 27.7 mag/arcsec 2 in the V-band, respectively. 

It is useful to compare this to the sample of late-type 
edge-on galaxies ob served by Dalcanton fe Bernstein! (200C) 
(their Fig. 3) and iDalcanton fc Bernstein (2002) (DB02, 
their Fig. 1), who typically probe up to ~ 5 mag below 
the central surface brightness of their (thin + thick) disks 
in the R-band. Their faintest contour would be located in 
between the first and second brightest contours shown in 
the right panels of Fig. [7] At least qualitatively, the surface 
brightness distribution of the remnants in our simulations 
resemble those observed by these authors. 

We further quantify the shapes of the isophotal contours 
of the rem nants by obtaining their photometry with the task 
ELLIPSE ([Jedrzeiew ski 1987) of the data reduction package 
This task draws an ellipse to approximately match 
an isophote and then expands the intensity alon g the e llipse 
as a Fourier series. According to iBender et all (|l988l ). the 
most significant non-zero component of this Fourier analysis 
is the a,4 parameter (corresponding to the cos(4#) term). 
Isophotes are then characterised as either disky (a,4 > 0) or 
boxy (d4 < 0). 

To mimic the observations, we have created artificial im- 
ages out of the simulated thick disks from an edge-on point of 



view, by binning a central area (~15x 15 scale- lengths of the 
initial primary disk) into 1024x1024 pixels. When running 
ELLIPSE on these images, we have allowed the geometric 
centre, ellipticity and position angle of the isophotes to vary 
freely, taking linear steps of 5 pixels along the semi-major 
axis. We have also made sure that our results are robust to 
the initial guesses for the values of the various parameters 
required by ELLIPSE. 

The left panel of Fig. [5] shows the 04 parameter as a 
function of isophotal surface brightness in the V-band for 
our "z=l" experiments. To avoid cluttering only experi- 
ments with spherical satellites on prograde orbits are ex- 
plicitly shown. The rest fall in the region delimited by the 
dashed envelopes. This figure shows that the isophotes go 
from more disky at higher surface brightness (i.e. in the 
central regions) to more boxy at lower surface brightness, 
presenting a mild trend with initial inclination of the satel- 
lite. 

|DB02| performed a similar analysis on their sample (see 
their Fig. 11). They find that inner isophotes with a sur- 
face brightness level of ~3 mag below the typical peak level 
for their sample (21 mag/arcsec 2 ) are disky, while outer 
isophotes (defined as those ~5 mag below the pe ak) are 
as likely to be boxy as disk}jf|. Similarly to |PB02| . in the 
right panel of Fig. [8] we plot the distribution of 0,4 (weighted 



2 IRAF is distributed by the National Optical Astronomy Obser- 
vatories, which are operated by the Association of Universities for 
Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 



J It is important to keep in mind that |PB02| use a different pro- 
cedure to quantify the isophotal shape, meaning that the values 
of their shape parameters are not directly comparable to ours. 
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Figure 8. Left: 124 vs surface brightness in the V-band for our 
"z=l" experiments of a spherical satellite on a prograde orbit. 
The dashed curves delimit the region where all experiments fall 
into. Note that the central surface brightness of our remnant disks 
are jJ.v.0 ~ 18 — 19 mag/arcsec 2 . Right: Luminosity weighted dis- 
tribution of the a4 parameter for all our "z=l" experiments. The 
dashed histogram corresponds to the £4 obtained by considering 
the region enclosed within an isophote with fj, < fig + 3, while the 
solid one to fig + 8 > fi > fio + 3. The arrows show the values of 
< 124 > for the coeval control simulation. 



by both errors and luminosity) for both inner and outer re- 
gions including all the "z=l" experiments. The inner region 
is defined to be within 3 mag from the peak surface bright- 
ness, as done by DB02. The outer region extends down to 
8 mag below the central surface brightness value. This fig- 
ure confirms that inner isophotes are disky while outer ones 
appear clearly boxy. This may su ggest t hat deeper photom- 
etry, beyond the limit reached bv lDBoj would be needed to 
detect the predominantly boxy shape of the contours in the 
outskirts of our remnants. 

The boxy nature of the outer isophotes, which is present 
in all our experiments, could in principle be used as a dis- 
criminant for the formation of thick disks via mergers such 
as those studied here. However, it should be borne in mind 
that the degree of boxiness in the remnants also depends on 
the initial structure of the host system. For example, studies 
which have a spherical centrally concentrated core compo- 
nent (bulge or cored dark matter profile) , p r oduce a remnant 
which is less boxy (Naab & Burkert 2003; iBournaud et al.l 
120051 ). This is because such sp herical components act in 
a sta bilising sense for the disk l|VW99l ; iKazantzidis et all 
120071 ). which therefore retains more closely its original mor- 
phology. This would imply that boxy isophotes are not nec- 
essarily direct evidence in support of the scenario proposed 
in this paper, but that they should be more prominent in 
the thick disks present in bulgeless galaxies. 
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Figure 9. Vertical luminosity profiles (integrating over all radii 
at each height) for each of our experiments. Fits using only one 
component (dotted lines) systematically underestimate the lumi- 
nosity far from the midplane. Two-component fits (solid lines) are 
clearly a better representation of the vertical structure of the our 
remnant disks. Note that, for clarity, profiles of inclinations of 30° 
and 60° include offsets of log(L) — 1 and log(L) — 2, respectively. 



3.4.2 Structural Properties 

We now describe in detail the structural properties of the 
remnant systems produced in our experiments. We first ad- 
dress their vertical structure and show explicitly the need 
to include two components (thin and thick). We then use 
the information obtained from the vertical decomposition to 
characterise the radial extension of both disk components. 
Finally, the spatial distribution of stars from both the pri- 
mary disk and the satellite are compared. 

Note that, as mentioned above, before measuring the 
disks, these have been properly centred and aligned, so the 



rotation axis defines the ^-direction. Furthermore, all the 
properties have been computed taking into account star par- 
ticles from both the host disk and the satellite. Their relative 
contribution has been appropriately weighted according to 
the initial M sa t, stars/ 'M host> disk ratio to account for the fact 
that the satellite particles have smaller masses. 

3.4.2.1 Vertical structure of the remnants Fig. 
shows the vertical luminosity distributions of the remnant 
systems for all our experiments. These have been obtained 
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Figure 10. Scale-heights of the final systems decomposed into 
a "thin" and a "thick" disk. Solid/dashed lines connect pro- 
grade/retrograde satellites. 
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Figure 11. Surface brightness profiles as a function of height at 
several projected radii for one of the "z=l" experiments. The thin 
remnant dominates the surface brightness in the central regions, 
while the thicker component is more dominant at the outskirts. 
Solid lines show the two-component fits at each radius. 



by including all stars at all radii within |z| < 3 kpc for 
"z=l", and within \z\ < 5 kpc for "z=0". The dotted lines 
correspond to a single-component sech 2 fit to the vertical 
profile. This figure shows that such a model clearly under- 
estimates the luminosity at large distances from the plane, 
highlighting the need for a two-component decomposition. 



L(R, z) = > L i (R) sech 2 



(5) 



where zo,i is a (luminosity weighted) exponential scale- 
height and Lo,i is the central luminosity (on the midplane) 
of each component (i — 1,2). As usual, fi{R,z) = 26.4 — 
2.5 \og[L(R, z)]. To compute the luminosity-weighted scale- 
heights we proceed as follows. First we fit independently the 
vertical brightness profiles in 1 kpc radial bins (within pro- 
jected radii R < 20 kpc for "z=0", and R < 10 kpc for 
"z=l") using two components. For each radial bin we allow 
the algorithm to find the best central luminosities and ex- 
ponential scale-heights using a Levenberg-Marquardt least- 
squares minimisation. The luminosity- weighted scale-height 
of each component is then the mean scale-height averaged 
over all radii and weighted by luminosity. 

The fits obtained in this way are shown in Fig. as 
solid curves. Clearly the vertical structure of our remnants is 
considerably better modeled by considering two disk compo- 
nents with different scale-heights and central surface bright- 
ness. In all cases, a thin disk is present after the merger with 
the satellite. 

Fig. [10] shows the luminosity- weighted scale-heights of 
each component of the remnant systems for all our experi- 
ments. Note that the thinner component has in all cases, a 
very similar (and only slightly larger) scale-height to that of 
the initial host disk. The scale-height of thicker component is 
clearly larger for encounters with higher orbital inclinations. 
This is because there is a significantly larger vertical kinetic 
energy associated to the satellites' orbital motion transferred 
to the disk. Spherical and disky satellite do not induce very 
different vertical heating on the disks. Note that less massive 
satellites produce final thick disks that are thinner. 

Fig. [TTJ shows surface brightness profiles as a function 
of height at several projected radii for one of the experi- 
ments. Two-component fits using the luminosity-weighted 
scale-heights described above are also included. This fig- 
ure shows that the remnant thin disk dominates the surface 
brightness at small radii. Note that at large radii (R = 9 
kpc, or jj, ~ fj,o + 6) there is an indication that the thick disk 
is flared, and no longer follows an exponential distribution 
with a constant scale-height at all radii. Such flar ed disk s 
have already been observed in previous studies (e.g. lQG86h . 
suggesting that flaring is a r ather generic characteri stic of 
disks heated by mergers fsee iKazantzidis et al.l [2007 . for a 
derivation of how the scale-height varies with radius due to 
minor mergers). 

3.4.2.2 Radial structure of the remnants Figures 
1121 and [131 show the mass surface densities of the simulated 
thick disks as a function of radius, for the "z=0" and "z=l" 
configurations, respectively. Each panel is for a different in- 
clination or different type of satellite (spherical, disky or 
half-mass) on prograde (red) and retrograde orbits (blue). 
The control disk galaxies are also shown at the initial (dot- 
ted curves) and final times (dashed curves) to calibrate the 
effect of the mergers against the intrinsic evolution of the 
host (which in all cases is negligible). In order to estimate 
the contribution of the satellite stars, the mass surface den- 
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Figure 12. Upper panels: surface density profiles of final thick disks as a function of cylindrical radius for "z=0", including stars within 
\z\ < 5 kpc. Red and blue lines correspond to prograde and retrograde orbits, respectively. The measurements include stars of both disk 
and satellite (solid) or only stars from the disk (dashed-dot). As reference, both initial (dotted) and final states (dashed) of the disk in the 
control model are also shown. Lower panels: Surface density profiles of regions dominated by thin (dashed; defined by \z\ < 0.5zo.thin) 
and thick (solid; Izq thin < I z I < 5 kpc) remnants, including the section used to compute the scale-lengths (heavy solid). 
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Figure 13. Same as Fig. 1121 but now for the "z=l" experiments, including stars within \z\ < 3 kpc in the upper panels. In this case 
the region dominated by thick remnants is defined as Izq thin < |z| < 3 kpc. The column "hm" corresponds to the lighter satellite (with 
half the mass). 
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periments. They are obtained by decomposing the final disks 
into thin and thick components. Solid/dashed lines connect pro- 
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sities of only host disks are plotted separately (dash-dotted 
curves) . 

For both "z=0" and "z=l" experiments, the surface 
density profiles show a mild dependence on orbital inclina- 
tion (upper panels), having slightly higher surface bright- 
ness in the outskirts for lower inclinations. This tendency 
is mostly due to the host disk material that is transported 
radially outwards during the merger by transfer of energy 
and angular momentum. 

As expected, the lighter spherical satellite (in the "z=l" 
experiment) produces a smaller variation in the surface den- 
sity at larger radius. In general the contribution of satellite 
particles to the spatial structure of final thick disks is very 
small. This is true at all radii except at the centre where 
the very dense spherical satellites accreted in the "z=0" ex- 
periment are not completely disrupted, retaining a core and 
giving rise to a small "bulge-like" component (see also Fig. 
[3]). The surface density profile does not show any strong de- 
pendence on orbit direction. 

The lower panels of Figs. [12] and [13] show the surface 
density profiles of regions dominated either by the thin or 
thick remnants (for prograde experiments only) according 
to the decomposition performed in £13.4.2.11 These regions 
are defined as \z\ < 0.5;zo,thm for the thinner component, 
and for the thicker one from lzo,thin upto 5 kpc for the 
"z=0" case (and upto 3 kpc for "z=l"). The scale- lengths 
of each component are computed by applying a linear fit to 
lnE(i?), avoiding non-axisymmetries associated to both the 
central regions and the very outskirts. Note that the linear 
fits consider a more extended region for thick remnants. This 
is to account for the dominance of thick remnants at larger 
radii (see Fig. Hip. 

Fig. [14] shows the final scale-lengths of the thick disks 



as a function of the initial orbital inclinations of the satel- 
lites. In all cases, the remnant thin disks have smaller scale- 
lengths than their thicker counterparts, and comparable to 
the initial values. Low inclination encounters induce larger 
thick-disk scale-lengths in comparison to higher inclinations. 
This is because in the former cases, the orbital energy of the 
satellite is deposited mostly into radial motions of the stars 
in the disk. Similar tr e nds a re observed for the galaxies in 
lYoachim fc Dalcantonl ((2006). However, this should be taken 
with great care because our comparison is to the remnant 
thin disk and not to the present-day thin disk of those galax- 
ies (because we do not model this). Furthermore we have not 
modeled the response of the remnant thin or thick disks to 
the formation of the new thin disk. 

We compute the total mass associated to each of disk 
component using the fits just derived. We find that the total 
mass associated to the remnant thin disk is ~15%-25% of 
the total stellar mass of the system for both "z=0" and 
"z=l" experiments. 

In general, the presence of a thin remnant after the 
merge r is in agreement with results by iKazantzidis et all 
(2007), although the total mass associated to this component 
is significantly smaller in our case (< 25% versus ~ 80%) . 
This is maybe due to the fact that IKazantzidis et all (|2007| ) 
do not follow the full merger event, but only let their satel- 
lites have one passage around their host disk, hence perhaps 
increasing its chances of remaining relatively cold. Note how- 
ever, that in their work this bombardment is repeated in a 
sequence using six different satellites. 

3.4.2.3 Distribution of satellite vs host disk stars 
Fig. 1151 shows the final relative number N S /(N S + Nd) of 
stellar satellite particles as a function of radius in the result- 
ing thick disks for mergers configured at "z=0" and "z=l". 
Recall that the number of satellite particles has been renor- 
malised according to N s = N aat ,stars X M sa t, s tars/M d i B k- 

As mentioned before, in the "z=0" experiments, spher- 
ical satellites have a higher mean density that the host disk. 
This causes the core to reach the host disk centre almost 
intact, representing ~50% of the total number of particles 
near the centre of the final thick disk. In comparison, this 
ratio drops to 5% for the disky satellites. At both intermedi- 
ate and larger radii the fraction of particles from the spher- 
ical satellites is roughly constant, independently of either 
orbital inclination or sense of rotation. On the other hand 
disky satellites are disrupted at large radii, where their de- 
bris is deposited. Furthermore, the lower the inclination, the 
higher the fraction at a given radius, as naturally expected. 

For "z=l" both spherical and disky satellites are com- 
pletely destroyed. The relative fraction of satellite stars in- 
creases with radius and, as expected, at the centre the rela- 
tive number of satellite particles is smaller for the lighter 
satellite. The observed trend with orbital inclination in 
spherical and disky satellites at "z=0" is confirmed for the 
"z=l" configuration. 

Fig. [16] shows the fraction of satellite particles plotted 
now as a function of the vertical direction at a radius 7?=4.5 
kpc, for spherical and disky satellites in the "z=l" experi- 
ments. The radius corresponds to 2.4 initial scale-lengths in 
this experiment. In this figure the distances from the plane 
are normalised by the scale-height obtained by fitting lo- 
cally a single sech 2 law to the vertical density distribution. 
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Figure 15. Relative number of satellite particles in the final 
systems as a function of cylindrical radius for the "z=0" (top) 
and "z=l" (bottom) experiments. 



Therefore Zo it is close to a luminosity weighted average of 
the scale-heights of the thin and thick disks given in Fig. 1101 
This figure shows that the fraction of accreted particles as a 
function of distance from the plane, when normalised by this 
scale- height, only depends on the mass ratio between the 
satellite and host. E.g., at z = 4Zo the fraction of particles 
reflects the mass ratio of the merger. The same behaviour is 
observed in the "z=0" experiments. 



3.4-3 Kinematical Properties 



The structural decomposition made in ij3.4.2. II should also 
be reflected in the kinematics of the stars in our systems in 
order to be physically meaningful. Fig. [17] shows this is in- 
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Figure 16. Relative number of satellite particles for our "z=l" 
experiments, as a function of distance from the disk plane. Dis- 
tances from the plane are normalised by the respective scale- 
height at a distance of R = 4.5 kpc (i.e. 2.4 scale-lengths of the 
original disk). Note that the scale- height here is obtained by fit- 
ting locally a single sech 2 law to the vertical density distribution, 
and hence it is close to a luminosity weighted average of the scale- 
heights of the thin and thick disk shown in Figure [Tul 



deed the case. Here we plot the z- velocity distribution within 
a spherical volume of 1 kpc radius centred at R ~ 4 kpc on 
the midplane for one of our experiments ("z—1", spherical 
satellite, prograde, 30°). The dashed curve shows that a sin- 
gle Gaussian (corresponding to a one-component system) 
misses the peak of the distribution highlighting the need 
for a second (colder) component. We therefore proceed to 
fit all velocity distributions (also the radial and azimuthal) 
with two Gaussians. We constrain the relative normalisa- 
tion of these by using the photometric decomposition from 
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Figure 17. Example of decomposition of the vertical velocity 
distribution into a cold component (associated to the remnant 
thin disk) and a hot one (the thick disk). The normalisation of 
each component has been fixed according to the photometric de- 
composition of £13.4.2.11 (see text). The contribution of the two 
components is separately shown. 



£ 13.4.2,11 which determines the relative number of stars from 
each component within a given volume. The solid curve in 
Fig. [T7]is an example of the quality of the fit obtained in this 
way, whose reduced-^ 2 (= 0.47) is lower than that obtained 
for a single Gaussian (= 0.68). 

The transformation of orbital energy into thermal en- 
ergy (random motions) in the disk is evident in Figs. I18l and 
1191 Here we show the radial, azimuthal and vertical velocity 
dispersions along with the mean rotational velocities of the 
thick disks present in our systems as a function of cylindrical 
radius. These quantities have been computed at each cylin- 
drical radius in concentric rings of 1 kpc width, including 
stars between \z\ < 3 kpc, for "z=0", and \z\ < 1.5 kpc, for 
"z=l" experiments. 

In all our experiments, the vertical and azimuthal ve- 
locity dispersions of the remnant thin disks are very similar 
to those of the initial host disk, and hence are not plotted 
for clarity. On the other hand, the radial velocity dispersions 
are generally larger by 5 — 10 km/s at all radii. 

Figures [18] and [19] show that the radial or and az- 
imuthal <70 velocity dispersions of the thicker component 
are larger for lower inclination orbits moving in the prograde 
sense. The opposite trend is observed for the vertical velocity 
dispersion az- This is as expected given our previous discus- 
sion on the evolution of the scale-heights and scale-lengths 
and their dependence on orbital inclination. Spherical and 
disky satellites give rise to similar velocity distributions. 

The resulting velocity dispersions or and az for retro- 
grade orbits are similar to the prograde cases. On the other 
hand, the azimuthal velocity distributions of the heated disk 
stars generally require an additional component to account 
for the contribution of the (accreted) counter-rotating stars 
(Villalobos & Helmi, in prep). The global velocity disper- 
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Figure 18. Final dynamical properties of the thick disk compo- 
nents for our prograde "z=0" experiments. These properties have 
been computed at t = 5 Gyr for stars within \z\ < 3 kpc in con- 
centric rings of 1 kpc width, using the decomposition described 
in £13.4.2. II As reference, the final state of the disk in the control 
model is also shown (dashed). 
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Figure 19. Final dynamical properties of thick disk components 
for our prograde "z=l" experiments. The column "hm" show the 
satellite with half the mass. As in the previous figure, the proper- 
ties have been computed within concentric rings of 1 kpc width, 
now for particles with \z\ < 1.5 kpc at t = 4 Gyr (except for the 
"hm" satellite, where t = 6 Gyr). The final state of the disk in 
the control model is shown by the dashed curves. 
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Figure 20. Mean azimuthal velocities of the "stars" in our experiments at several heights above the plane as a function of radius for 
our "z=l" experiments in prograde orbits. Thick lines include particles from both the heated disk and satellite, while the thin lines 
correspond to the disk in the control model (the highest z-range is not shown in this case because there are too few disk particles at that 
distance from the midplane). The same behaviour is observed in our "z=0" experiments (not shown). 



sion (that would be obtained by imposing a single Gaussian 
for the thick component) would be in this case significantly 
larger at large radii for retrograde orbits, except for those 
with high inclination. 

In general, the mean rotational velocities v$ of the thick 
disks differ noticeably from the coeval control simulation in 
all cases, by dropping ~60 km/s although with a mild de- 
pendence on the orbital inclination. Low inclination encoun- 
ters produce thick disks that rotate slower, implying larger 
asymmetric drifts. This is also evidenced by their larger ra- 
dial and azimuthal velocity dispersions, as discussed above. 
The mean rotational velocity of thick disks which are the 
result of encounters with satellites on counter-rotating or- 
bits is lower due to the contribution of the accreted stars, 
particularly at large radii. 

The mean rotational velocity also shows noticeable dif- 
ferences with inclination when it is measured away from the 
midplane of the thick disk. In Fig. [20] we plot v$ for the 
prograde experiments at different heights above the plane, 
without making a distinction between the thin and thick 
disk components. Note that for \z\ > 1 kpc we are really 
measuring the kinematics of the thick disk component since 
the contribution of the remnant thin disk is negligible. Satel- 
lites with lower initial orbital inclinations induce a rotational 
lag, whose magnitude increases with height above the plane. 
This is because such satellites are more efficient in heating 
the disk radially at every height, leading to a larger asym- 
metric drift. 
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Figure 21. Final az I^R of thick remnants as a function of dis- 
tance on the plane for different initial stellar distributions and or- 
bital inclinations of the satellites. The control disks have a steady 
uz/cr ~ 0.7. Only prograde orbits are shown. 



In Fig. [2T] we plot the ratio az I or of the thick disk 
component as a function of radius for different prograde ex- 
periments. Recall that the initial (and the control) disk has 
a (nearly) constant cfz/&r ~ 0.7. This figure clearly shows 
that az/o-R can be used as a discriminant of the initial incli- 
nation of the satellite. The reason for this strong dependence 
on inclination is essentially due the fact that a satellite on 
a highly inclined orbit will induce a much larger change in 
az at large radii than one on a co-planar orbit. 



4 DISCUSSION 

Previous works investigating the response of a host disk 
galaxy to one or more infalling satellites provide us with a 
valuable descrip t ion of the dynam i cal effects invo l ved in the 
process (IQG86I; iToth fc Ostrikeil 1 19921; IQHF91 ; IWMH9& 
lHC97l ; IVW99l ; lBenson et al.ll2004h . 1110971 find that a massive 
satellite (30% of Mdisk)t as ^ decays due to dynamical fric- 
tion, can tilt the orientation of the disk up to 10° and cause 
warping. These effects illustrate the strong transfer of angu- 
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lar mo mentum from the infalling satellite to the host disk. 
IVW99I note that satellites in prograde orbits mostly increase 
the disk heating. On the other hand retrogr ade orbits are 
more efficient in tilting the disk orientation. iBenson et al.l 
show that more massive and more concentrated satel- 
lites increase the difference between the amou nt of d isk heat- 
ing caused by prograde and retrograde orbits. QG86 indicate 
that most of the kinetic energy of the infalling satellite is de- 
posited on the plane of the disk and a small quantity heats 
up the vertical motions of disk particles. This shows that ve- 
locity dispersions in t he host d isk are not increased isotrop- 
ically by the satellite. IQHF93I note that most of the heating 
on the plane of the disk is caused by spiral arms stimulated 
by the decaying satellite. These arms also transfer angu- 
lar momentu m radially outwards, expanding the disk. lQG86l 
and QHF93 also show that the vertical structure of the disk 
is not uniform across the disk. Instead, the scale-height of 
the host disks increases at larger radii. Our simulations are 
able to confirm most of the aforementioned dynamical ef- 
fects but also show a significantly larger tilting than pre- 
viously found, both for the prograde and retrograde cases. 
This can be traced back to the fact that in our experiments, 
the accreted satellite is initially much more massive and is 
launched from a significantly larger distance (see H3.3|l . 



4.1 Choice of massive satellites to heat disks up 

Cosmological simulations show that massive mergers like 
those in our "z=l" experiments are likely to have happened 
during t he life-time of a Mi l ky W ay-sized host system. For 
instance, iKazantzidis et all (|2007T > estimate the number of 
massive subhalos accreted since z ~ 1 as at least one object 
with a mass ~Mdisk an d five objects more massive than 20% 
Mdis k- This is consistent with simulations by other groups 
(e.g. IStoehrll200d |Pe Lucia fc HelmilbOOsl). and supported 
by the detailed study by Stewart et al.l ( 20071 ). Even though 
massive mergers may be less frequent they are able to reach 
the centre of the host system thanks to dynamical friction, 
causing important changes in the structure and kinematics 
of the host disk. In principle, this could imply that a Milky 
Way-like disk would experience only one severe change of 
orientation since z ~ 1 and also that only one massive satel- 
lite would be needed to heat up a pre-existing disk to give 
rise to a thick disk. 



4.2 Observations of thick disks 



iPohlen et al.l l|2004 ) have carried out photometric thick/thin 
disk decompositions and characterised properties such as 
scale-length and scale-height for a sample of eight SO 
edge-on galaxies by fitting 3D disk models. The authors 
find that the mean scale-height of the thick disk is be- 
tween 2.4 and 5.3 times larger than that of the thin 
disk; while the mean scale-length is about twice, rang- 
ing from 1.6 to 2.6. In general these values are consis- 
tent the results obtained by other authors for SO , Sab 
Sb, Sbc, Scd and Sd galaxies (e.g.. Ivan der KruitJ 



de Griis fc van der KruitJ Il996l; Ide Griis fc Peletieri 



1984 



1997 



Pohlen et al.ll2004l; lYpachim fc Dalcantonl 120061), includ ing 
the Milky Way (|Oihall200ll ; lLarsen fc Humphrevsll2003l ). // 

we assume that the structure of our simulated disks are not 



strongly affected by the formation of a new thin disk, and 
that these new thin disks will follow the same distribution 
as the remnant thin disks in our models, we can take the 
ratio of scale-heights obtained for our experiments at the 
final time at face value. Typically there is an increase in 
the scale-height by a factor 3-6, while the scale-lengths 
are only slightly larger (with final-to-initial ratios between 
1 and 1.4). These values appear to be in agreement with the 
range of ratios observed in general in spiral galaxies. How- 
ever, this should be taken with great care because of the 
strong assumptions just made. 

Most of our thick disks are significantly flared in the 
oute r regions. This has also been found in o ther studies 
(e.g. IKazantzidis et all [20071 ; iRead et all 120081 ). with vary- 
ing strengths, depending on the initial configuration of the 
system. For example, galaxies with a bulge or embedded in 
a cored dar k-matter halo are more stable, and suffer from 
less flaring dNaab fc BurkertJ 120031 ; iBournaud et all 120051 ; 
IKazantzidis et al.ll2007 ). This implies that we might expect 
to find in nature thick disks with varying degrees of flar- 
ing. However, quantifying this requires reaching extremely 
low surface brightness levels, approximately 6 — 7 magni- 
tudes below the peak brightness of the thick disk. This is 
very challenging and has not yet been achieved in studies 
of surface photometry, which typically reach about 5 mag- 
nitudes below the peak value of the thin + thick disks. For 
example, Neeser et al (2002) detected a thick disk in a low- 
surface brightness galaxy. They reach a surface brightness 
level in the R-band of 29 mag/arcsec 2 . This is almost 7 mag- 
nitudes below the central thin + thick disk value, and a very 
marginal indication of a flare is apparent on the east side of 
the thick disk at a distance of ~ 3 scale-lengths. Morrison et 
al. (1997) in their study of NGC891 (a late type disk with 
a small bulge) do not find evidence for flaring in the thick 
disk, but "only" reach down to 26 mag/arcsec 2 in R, which 
is about 6 magnitudes below the central surface brightness 
of the dominant disk component. 

iDalcanton fc Bernstein! (|2002h have found a vertical 
colour-gradient in their sample of galaxies, in the sense that 
the thick disks are redder, and consistent with a relatively 
old (> 6 Gyr) stellar population. In the scenario we envisage 
for the formation of thick disks we would expect two sources 
of color-gradients. The first would be due to the transition 
from the population of young thin disk stars to the thick disk 
stars formed in situ, i.e. originating from a pre-existing thin 
disk. Unfortunately, because we do not model the formation 
of a new disk from freshly accreted gas, we cannot make 
a quantitative comparison. The second source of a change 
in stellar populations would be due to the difference be- 
tween accreted satellite stars and those formed in situ in the 
(heated) disk. Our Figures [Pol and [TBI show that the fraction 
of satellite particles is generally quite low. Only at distances 
above the plane of approximately 6 thick disk scale-heights, 
the number of accreted stars is comparable to that of the 
heated disk. This implies that this second gradient should 
be weaker, and is unlikely to be detected in current observa- 
tions because at these heights, the typical surface brightness 
levels are exceedingly low (see Figure ITT]) . 

Kinematically, the values of the velocity ellipsoids 
of our thick disks, measured at R ~ 2ARd, are in good 
agreement with the one observed at the solar radius in 
the Milky Way (o\R,o-,/,,erz)~(65,54,38) km/s (e.g., see 
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Lavden et all Il996l; Ichiba fc Beersl l200ll; ISoubiran et all 



20031; lAlcobe fe Cubarsil 120051 ; Ivallenari et all l200rf 
Veltz et al.ll2008h . The values of az/o-R measured by most 
of these authors are typically ~ 0. suggesting that the 
thick disk of the Galaxy could have been produced by 
a merger of intermediate inclination (see Fig. 1 2 1 fl . The 
differences between the mean rotational velocity of the final 
thick disk and that of the remnant thin disk are A^~40-50 
km/s. Note that these values cannot be interpreted directly 
as the observed rotational lag between the thin and thick 
disks in the Milky Way since we do not include gas in our 
simulations that later on can collapse and actually form a 
new thin disk (see ^4.4. 1 p . 

4.3 Our simulations in the context of M31 

It is important to notice that in our experiments config- 
ured at "z=0" the initial properties of the host disk, before 
the merger with the satellite, resemble the properties of the 
current thin disk of the Galaxy. This implies that the final 
thick disks generated in our experiments are too massive 
and do not represent direct analogues of the thick disk of 
the Milky Way. However this set of simulations may be use- 
ful for understanding the evolution of M31, which is likely 
to have experienced recent merger events as may be inferred 
from the rich and complex structures present in its outskirts 
jlbata et al.l 120071. and references th erein) . 

For example, Ibata et al.l l|2005l ) discovered an extended 
and clumpy disk- like structure around M31 with a scale- 
length similar to that of the main disk, rotating ~40 km/s 
slower and with a rather low velocity dispersion of ~30 
km/s. Its stellar population has homogeneous kinematics 
and abundances over the entire region where it is observed, 
which suggests that it was formed in a single global event. 
However, it is not straightforward to link the extended disk 
of M31 to the thick disk modeled here, mainly because of 
the much higher velocity dispersions of the final thick disks 
in our simulations. 

The giant stream of M31, thought to originate in the 
disruption of a satellite with mass ~10 9 M^, and t he so- 
called eastern and western "shelfs" (|lbata et al.l l2007), may 
be related given the similarities of their stellar populations. 
In the context of our simulations, this would be quite natu- 
ral. Analogous structures are observed as the satellite is dis- 
rupted as shown in Figs. [5] and E] for prograde orbits in face- 
on vi ew at 1.5 Gyr (see also iFardal et al.ll2007l ; iMori fc Rich! 
2008). It is interesting to note that much sharper and longer 
lasting shells are generated by disky satellites compared to 
spherical ones. 



nor star formation which may affect the interactions and 
the final thick disk's properties. This is potentially the most 
crucial simplifying assumption in this study since disk galax- 
ies were presumably m uch more gas rich at high redshift 
l|Robertson et al.ll2006l ). 

The lack of gas physics implies that the modeled disks 
do not grow in stellar mass or in size during the times- 
pan of the merger (except of course through the dynami- 
cal heating processes described above). Furthermore, it is 
likely that such mergers would trigger a (strong) burst of 
star formation. These new stars would be relatively old at 
the present day and located in a thinner structure. On the 
other hand, some orbital energy deposited by the infalling 
satellite into the gas could be radiated away from t he system , 
reducing the dynamical da mage done to the disk (|QHF93h . 
For example, in rece nt work iHopkins et all (|2008l ) based on 
lYounger et alj (2008) suggested that the change in the struc- 
tural parameters depends on the fraction of gas f g available 
as 8H ~ (1 — f g )5H,, where 5H* corresponds to the scale 
change in the purely dissipationless case. 

Any remaining (and presumably heated) gas would 
eventually cool and settle down to form a new thin disk. 
This slow accumulation of gas on the midplane should 
induce a contraction of the thick disk. For example, 
lElmegreen fc ElmegreerJ (|2006l ) estimate that this contrac- 
tion leads to a decrease in the scale-height of the contracted 
thick by ~40% and an increase of the velocity dispersion 
by ~50%. Furthermore, the accretion of fresh gas from the 
inter-galactic medium is likely to also be important, and 
will lead to further changes to the properties of the merger 
products studied here. 



4-4-2 Time- dependence of gravitational potential 

In general, the structure of a dark matter halo evolves with 
time through mergers and slow accretion. However, in our 
simulations we have neglected any cosmological evolution of 
the structure of the host halo during and after the merger 
with the satellite. This simplification may be justified by re- 
cent studies fe.g.. IWechsler et alj |2002; Ro mano-Diaz et al.l 
2006), which have shown that the structure of dark halos 
is very stable within the scale radius r s after the phase of 
active mergers, which for a disk galaxy must have ended at 
redshifts ~ 0.5 — 1. This is indeed the region that we fol- 
low dynamically in our simulations, after the initial decay 
of the satellite due to dynamical friction, which lasts typ- 
ically less than 1 Gyr. Therefore, the final thick disks are 
well within this scale radius, being r s ~ 6 and ~ 11 times 
the final scale- lengths in the "z=0" and "z=l" experiments, 
respectively. 



4.4 Caveats 

4-4-1 Lack of gas physics in our simulations 

In this paper we have focused on the collisionless interactions 
between a disk galaxy and a satellite, with both dark mat- 
ter and stellar components, without including gas physics 

4 Although Veltz et al find uz/ur ~ 0.9, which could be due 
to the fact that the assumption of isothermality for the velocity 
distribution is not valid. 



5 SUMMARY AND CONCLUSIONS 

We have performed numerical simulations of the heating of 
a disk galaxy by a single relatively massive merger. These 
mergers lead to the formation of thick disks whose charac- 
teristics are similar, both in morphology as in kinematics, to 
those observed in the Milky Way and other spiral galaxies. 

The simulations explore several configurations of the 
progenitor systems whose properties have been scaled at two 
different redshifts in order to study the formation of thick 
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disks at different epochs. The satellites have total masses of 
10% and 20% that of the host galaxy and have been modeled 
self-consistently as a stellar component immersed in a dark 
matter halo. The stellar components have either a spherical 
or disky distribution. The satellites have been released far 
away from the host disk, with initial orbital parameters that 
are consistent with cosmological studies. Additionally, three 
different initial orbital inclinations of the satellites have been 
studied in both prograde and retrograde directions with re- 
spect to the rotation of the host disk. 

We find that as the satellite galaxies spiral in through 
dynamical friction, significant asymmetries are visible, both 
in the host disk and in the satellite debris. Particularly inter- 
esting are the low surface brightness shells, especially visible 
in the outskirts of the final thick disks, that last for about 1.5 
to 2 Gyr after the merger has been completed. These shells 
acquire relevance in the case of Andromeda where according 
to recent studies a couple of these features are likely asso- 
ciated to the even t that also gave rise to the giant stream 
|lbata et al.ll2007l and references therein). 

Despite the relatively large mass ratios, the infalling 
satellites do not fully destroy the host disk, but merely heat 
it and tilt it. The host disks are found to change their ori- 
entation both for the prograde and retrograde encounters. 
Furthermore, a remnant thin component containing between 
15 and 25 per cent of the total stellar mass of the system 
is present at the final time in all our experiments. This pre- 
diction of the minor merger model might be potentially rel- 
evant to understand the presence of a very old thin disk in 
the Milky Way. 

The scale-lengths of the final thick disks are slightly 
more extended than those of the original host disk while 
the scale-heights are between three and six times larger, de- 
pending on the initial inclination of the satellite. The scale- 
heights have also increased in proportion to the inclination 
of the encounter, and the outer disks are noticeably flared. 
If this is the case for the thick disk of the Milky Way, part of 
the flared material could be (confused with) the Monoceros 
ring. 

In our simulations, the outer isophotes of the final thick 
disks (measured at surface brightness levels >6 mag below 
the central value) are consistently more boxy than the inner 
ones. The eventual detection of such degree of boxiness, es- 
pecially for bulgeless galaxies, would provide support for a 
formation process as that modeled here. 

Interestingly, the fraction of satellite particles at a given 
galactic radius as a function of height above the plane, when 
normalised by the (luminosity-weighted) scale-height, only 
depends on the mass ratio between the satellite and host and 
not on stellar morphology of the satellite or type of orbit. 
For instance, at a distance of 4 scale-heights the fraction of 
satellite particles reflects the mass ratio of the merger. 

We find that satellite stars do not dominate the lumi- 
nosity of the thick disk until rather far above the midplane. 
In this sens e, the existence of a counter- rotating thick disk, 
detected bv lYoachim fc Dalcantonl (|2005h only ~2 thick-disk 
scale-heights above the midplane, can only be explained in 
the context of our models if the (young) thin disk formed 
from freshly accreted counter-rotating gas. The remaining 
possibility is, of course, that the thick disk formed exclu- 
sively by direct accretion of stars from an infalling satellite. 
Relatively fast rotating thick disks (like the one of the Milky 



Way) may be more easily explained by disk heating forma- 
tion, since a random distribution of accreted satellites would 
seem to have less chance of producing thick disks with strong 
coherent rotation. 

// taken at face value the velocity ellipsoids of the simu- 
lated thick disks are in good agreement with observations of 
the Galactic thick disk at the solar radius. The rotational lag 
may also be consistent with observations. These statements 
are however only valid if we neglect further significant evo- 
lution due to the formation of a thin disk component from 
freshly accreted gas. The observed trend of the ratio oz I or 
with radius in the final thick disks is found to be a very 
good discriminant of the initial inclination of the decaying 
satellite. In the case of the Milky Way, the observed <jzIo~r 
at the position of the Su n is ~ 0.6 (e.g., IChiba fc Beersl 
l200ll : IVallenari et al.ll2006l ). suggesting that the thick disk 
of the Galaxy could have been produced by a merger of 
intermediate inclination. Measurements of the mean rota- 
tional velocity in the final thick disks, at several heights from 
the midplane, indicate that satellites with lower initial incli- 
nations are more efficient in introducing asymmetric drifts 
dependent on height. This implies that the possible exis- 
tence of vertical gradients in th e mean rotational velocity 
in the thick disk of the Galaxy (jGirard et al.ll2006l ) would 
also favour mergers with either low or intermediate orbital 
inclination. We defer to paper II a more detailed analysis of 
the phase-space structure of the merger product. We expect 
this will lead to new constrains on the mechanism described 
here for the formation of the Galactic thick disk. 
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APPENDIX A: SETTING UP THE INITIAL 
CONDITIONS FOR THE HOST AND 
SATELLITE SYSTEMS 

Al Main disk galaxy 

The main disk galaxy is a self-consistent two-component sys- 
tem, containing a dark matter halo and a stellar disk. 



A 1.1 Dark Matter Halo 

The dark matter (DM) ha los in our simulations follow a 
NFW mass density profile l|Navarro. Frenk. fc White! 1 19971 . 
hereafter NFW): 



pNFw{r) = 



(r/r 3 )(l + r/r s ) 2 



(Al) 



where p s is a characteristic scale density and r 3 a scale ra- 
dius. The advantage of using this density profile is that it 
is consistent with cosmological simulations, and its evolu- 
tion with redshift (o r time) is relatively well-known (e.g., 
IWechsler et"aT]|2002l ). This implies that it is easy to re-scale 
its properties to study the formation of thick disks at red- 
shifts greater than zero. 

In this paper we adopt a flat cosmology defined by 
fi m (z = 0) = 0.3 and JIa = 0.7 with a Hubble constant 
of H{z = 0) = 70 km/s/Mpc. 

The virial radius of the halo R v i r (z) is defined as the 
radius within which the mean density is A„i r (z) times the 
critical density p c (z) of the universe at a given redshift: 



Am 3 
M vir (z) = — A vir (z)p c (z)R vir (A2) 

where the virial overdensity A„i r (z) is taken from the solu- 
tion to the dissipationless collapse in the spherical top-hat 
model. Its value is 187r 2 for a critical universe but has a 
dependency on cosmology. In the case of flat cosmologies, 
A vir (z) « (187r 2 + 82x + 39a; 2 ), where x = Q(z) - 1, and 
f2(z) is defined as the ratio between mean matter density 
and critical density at redshift z. Another important related 
quantity is the c oncen tration c defined as c = R V i r /r 3 . From 
IWechsler et ail (|2002ft we take the relation linking M V i r to 
the concentration parameter c at redshift z=0 as: 

-K^)" 013 - (A3) 

We follow IWechsler etaD l|2002h to scale both the virial 
mass of the halo and its concentration as a function of red- 
shift: 



M V ir(z) = M vir (z = 0) exp(— 2a c z), 



c(z) 



c(z = 0) 
1 + z 



(A4) 
(A5) 



where a c is a constant defined as the formation epoch of the 
halo, taken as a c = 0.34. In practice this means that the 
structure of the halo of the main galaxy at any redshift is 
fully determined by imposing only a value for the virial mass 
at redshift z=0. The values of the halo parameters used in 
our simulations are included in the Table [1] 

Since the mass of a NFW profile formally diverges with 
radius, we introduce an exponenti al truncation starting a t 
Rvir and decaying on a scale rdec l|Springel fc Whitelll999h : 

= <i+~c7 o-y^ (-^) ^ > ^ (Ae) 

where rdec is a free parameter. By requiring continuity at 
Rvir between Eqs. (|A1|I and (|A6|I . and also between their 
logarithmic slopes, the exponent e is computed as: 

£ = _1-3C + Rvir (A7) 

1 + c r dec 

Note that for rdec = 0.1R V i r the total mass of the halo 
becomes ~ 10% larger than M V ir. We define the maximum 
extension of the halo as R m ax = Rvir + 3rd ec . 

We also allow the contraction of the halo in response 
to the formation of a stellar disk in its central part 
l|Blumenthal et al.lll986l : iMo et al.lfl99§ ) The adiabatic con- 
traction first assumes that the gas (that later forms the 
disk/bulge) is distributed in the same way as the dark mat- 
ter. Then both the spherical symmetry of the halo and also 
the angular momentum of each dark matter orbit are con- 
served during the contraction, i.e.,: 



nMiin) = r f M f (r f ) 



(A8) 



where r» and r/ are the initial and final radius of a shell of 
dark matter, Mi is the initial total mass (distributed accord- 
ing to an NFW profile) and Mf is the mass distribution after 
the disk has been formed, and also includes the contribution 
of the disk. Therefore, 



M f (r f ) = M d (r f ) + (1 - m d )Mi(n) 



(A9) 



where m d = Mdisk/Mhaio- The final dark matter distribu- 
tion of the adiabatically concentrated halo will be: 
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M halo (r) = M f {r) - M disk {r). 



(AlO) 



Now that the mass distribution of the halo component 
has been defined, it is straightforward to initialise the posi- 
tions of the particles in our halos. As the next step, the veloc- 
ity of each particle is computed from the distribution func- 
tion (DF) associated to the adia batically contracted mass 
density profile phaio(r). We follow IKazantzidis et al.l (|2004l ) 
and compute numerically the DF that, in general, is given 



f(Q) = 



(fphalo dip 



dip 2 VQ^P 



+ 



dphal. 

VQ V # 



i/i=0 



(All) 



jBinnev fc Tremaindll987l ) where Q = ip-v 2 /2; ip = -$(r) 
is the effective gravitational potential (including the disk); 
and v is the vel ocity of each parti cle. Finally, we use the re- 
jection method (|Press et al.ll 19921 ) to generate the velocities 
for our particles. 

A 1.2 Stellar Disk 

The disk co mponent is c onstructed f ollowi ng the procedure 
outlined by IQHF93I and iHernquistJ (|l993l ) which consists, 
briefly, in initialising particle positions according to a den- 
sity profile of the form: 



pd(R,z) = 



M d 



■ ,;; ,' ( ' X|> 1 /»'/.) S6Ch ( 2s,; ) 



(A12) 



where M d is the disk mass, Rd is the exponential scale- 
length, and zq is the exponential scale-height. 

The velocity components vr, Vj, and v z of the disk par- 
ticles are calculated from moment equations of the colli- 
sionless Boltzmann equ ation (CBE) supplemente d by ob- 
servational constraints (IBinnev fc Trem ainc 1983). We as- 
sume that locally (at each point in the disk) the velocity 
distribution can be approximated by a Maxwellian, whose 
parameters are set up as follows: 

• The radial velocity dispersion v 2 R (R) oc exp(— i?/i?n). 
Thi s is motivated by observat i ons of external galax - 
ies (Ivan der Kruit fc Searlelll981al ; iLewis fc Freemanll 19891) . 
The normalisation constant is set by adoptin g a certain value 
of the stability Q-parameter l|Toomrd [l964) at a particular 
location in the disk. In this paper Q = 2 at R = 2ARd, 
which for a Milky Way-like disk corresponds to the solar 
radius. 

• The vertical velocity dispersion v 2 (R) — 2nGE(R)zo, 
following the isothermal sheet model. 

• The dispersion in the azimuthal direction is obtained 
by using the epicyclic approximation (Bin nev fc Tremaind 
1 19871 ) al(R) = v^(R)k 2 (R)/4Q, 2 (R), where k and Q, are the 
epicyclic and angular frequencies, respectively. The mean 
values of the azimuthal Gaussian distributions are calculated 
from the second moment of the CBE, 



vf{R) = v 2 R {R) 



\R) 



R 



4fi 2 (i?) R D 



■v 2 (R), 



where v c (R) — R f2(i?) is the circular velocity considering 
all the components of the system. 

Note that velocities derived from the CBE are close but 
not identical to the ones derived from the DF of the disk. 
Unfortunately, the DF is unknown for the disk in Eq. <| A12[) . 
Therefore, we can expect some initial evolution in the disk 
properties. As shown in ^ A31 this evolution is indeed mini- 
mal. 



A2 Numerical Parameters 



For a spherical non-rotating system. 



The A ^-body systems are evolved using Gadget-2.0 (|Springell 
a well documented massively parallel TreeSPH code. 
Depending on the system under study, this code has to be 
provided with suitable values for the so-called numerical pa- 
rameters, being these: the number of particles N to represent 
a given component in the system; the softening e of gravita- 
tional forces to avoid strong artificial accelerations between 
particles passing close to each other; and finally the timestep 
At, that controls the frequency at which positions and veloc- 
ities are computed for each particle. In general, these three 
parameters set the mass, spatial and time resolution in a nu- 
merical simulation. At the moment of defining N, e and At 
the usual problem is that they are interrelated in a compli- 
cated way. For instance, N will depend on the available CPU 
power to run the simulations; e will depend on both A^ and 
the mass distribution of the system to be simulated; and At 
will depend on the smallest spatial resolution that is possible 
to resolve, e, and again on the available CPU power. The op- 
timal choice of these parameters will establish a compromise 
between quality and efficiency in a numerical simulation. 



A2.1 Number of particles 

Tables [T] and [2] list the numbers of particle s used fo r each 
component in our simulations. As shown by WMH96, using 
self-consistent simulations of an isolated halo-disk system, 
large numbers of particles in the halo are needed to sup- 
press the formation of bar perturbations in the disk. This is 
because large N^aio decreases t he grainin ess of the poten- 
tial, which bars are seeded from. IWMH96I suggest the use of 
~ 500000 particles in the halo in order to smooth out the 
potential for time scales comparable to the orbital decay of 
satellites in our simulations. 

For our purposes, bars are an unwanted effect because 
they represent an additional source of disk heating, besides 
the one of interest here. Although the complete elimination 
of bar formation in a self-consistent simulation is difficult, 
its effect on the disk can be constrained by evolving the 
main disk galaxy in isolation, for the same timescale as the 
merger simulation. 

The number of particles in the host disk Nhost,disk 
10, and is similar to previous studies on disk heating by 
mergers with satellites. The satellites are modeled with a rel- 
atively large number of particles (particularly in comparison 
to previous works) to study the distribution of the debris, 
which is the focus of Paper II. In all cases we can follow 
accurately the structure and evolution of each component 
during the simulations. 



Simulations of minor mergers. I. 23 



o.i 



10 2 



Halos 


Disks and Rulges 


_ Qhost @> 7, = 


_ O d host • z =° 


Ohost 8 z=l 


0D host z=l 


Aspher sat @ z = 


AB sphcr sat @ z = 


" Aspher sat @ z=l 


:AB spher sat @ z=l 


_ □ disky sat @ z=0 _ 


Z □ D disky sat z = 


K 5'4 disky sat 'S 1 z 1 . 


^ D disky sat @ z ::: 1 






% : 




1 O 


n 



0.02 0.04 0.02 0.04 

r 6/ r i/2 

Figure Al. Numerical softening as a function of the mean dis- 
tance from each particle to the 6th closest neighbour in units of 
the half-mass radius fi/2' Straight lines ar e extrapolations to- 
wards smaller re taken from the results of lAt hanassoul a et al.l 
[see Fig. 11 in that work]. Dashed, solid and dot-dashed 
lines show the values for a homogeneous, Plummer and Dehnen 
sphere (with 7 = 0), respectively. The symbols indicate the values 
of softenings explored in this paper for each component. Darker 
symbols show the optimal softenings that produced the best sta- 
bility in each case. 



A 2. 2 Softening 

Many studies have been carried out on how to choose the 
optimal gravitational for ce softening e in order to faithfully 
represent the system (see lSellw ood 1987 1 for an excellent re- 
view) . We use here the prescription by lAthanassoula et al.l 
(2000). These authors present a simple method to estimate 
e for arbitrary mass distributions as a function of the num- 
ber of particles. The optimal softening e opt is the one that 
minimises the error in the forces between particles in a sys- 
tem with a given mass distribution. They find a correlation 
between e op t and the distance re,mean from every particle to 
its sixth closest neighbour, which is defined as: 

ni.mean = ^ 7^ (A13) 

where re,mean depends on both the number of particles 
and the mass distribution. Fig. IA1I shows func- 
tion of r§, mean for three mass distributions discussed by 
lAthanassoula et alj (|2000l ) in order of increasing density: 
homogeneous sphere, Plummer profile and Dehnen model 
(with 7 = 0). To estimate the e op t for the systems of our 
simulations the procedure followed is: 1) compute re, mean 
for each of our components; 2) compare the central density 
of our components to the ones of the homogeneous, Plum- 
mer and Dehnen spheres. By doing so the optimal softening 



for each component can be constrained within a range on 
the plane re — e opt ; and 3) e opt is found by running a few 
simulations with a set of e within this range and choosing 
the one that offers the best stability. Specifically, a soften- 
ing is considered optimal when each component of both the 
host and satellite systems present the least evolution in their 
structure and kinematics (in the case of the host minimising 
the effect of non-axisymmetries) during the amount of time 
required for the satellite to sink and become fully disrupted. 

In Fig. IA1I darker symbols show the adopted values of 
t op t for each component at every redshift. Note that for each 
component the values of e op t are well constrained on the 
re — e op t plane, facilitating the extrapolation of these op- 
timal values to similar systems with a different number of 
particles. For each component, the adopted values of e op t are 
listed in Tableland Tabled 

The distinctive location of the optimal softenings on 
the re — e op t diagram basically depends on both the central 
concentration of the systems and on the number of particles 
used to model them. For instance, halos of hosts and satel- 
lites at "z=0" and "z=l" always lie below the Dehnen model 
because they are more centrally concentrated (even more so 
when the adiabatic contraction is taken into account). On 
the other hand, the separation along the r6-coordinate be- 
tween the optimal softening of the halos of hosts and satel- 
lites reflects the difference in the number of particles used 
to model them. In this sense it is easy to associate systems 
with a larger number of particles to smaller mean distances 
between particles and viceversa, given that the systems are 
compared in a normalised scale. 

Also note that the disky satellite requires a halo that 
is better resolved at "z=l" than for a spherical satellite in 
order to reach the best stability. This can be explained by 
the fact that a better resolved centrally concentrated re- 
gion of the halo is ab le to inhibit the formati on of non- 
axisymmetries (e.g.. see lAthanassoula et alJ feoOS). 

In order to check the robustness of our choices of both 
number of particles and softenings, we have followed the sug- 
gestion by the referee and also simulated one of the "z=l" 
experiments adopting a the same mass and softening for the 
stars in the satellite and in the host disk stars. Reassuringly 
we found practically no difference in the global properties of 
the final thick disk in comparison to our "standard" choice 
of numerical parameters. 

A 2. 3 Timestep 

The timestep At has been defined for our simulations ac- 
cording to the standard criterion of Gadget-2.0. This means 
th at the timestep for each particle is calculated as At = 
y^2ne/\a\, where n is a dimensionless parameter controlling 
the accuracy of the timestep criterion, e is the softening 
associated to each particle, and a is the gravitational ac- 
celeration suffered by each particle. At is also limited by 
a maximum value in order to prevent particles having too 
large timesteps. The maximum timestep is defined as a few 
percent of the timescale t c = 2ne/V c (e) calculated for the 
component that has the smallest e in the system, where V c 
is the circular velocity. This ensures us that we follow the 
evolution of even the smallest components in the system 
with enough time resolution. We have set r\ = 0.025 and the 
maximum timestep to 0.25 Myr, which give us typical con- 
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in isolation only after 9 Gyr and 7 Gyr for the configurations 
at "z=0" and "z=l", respectively. 

Fig. IA2l shows how the initial properties of the host disk 
at "z=0" change after being relaxed within a fixed halo for 
1 Gyr, and after 5 Gyr in the live halo. Its properties are 
measured in concentric rings of 1 kpc of width, including 
particles out to ~15 initial scale-heights. The surface den- 
sity profiles E(-R) (top left panel) indicate that the scale- 
lengths of the disks (given by the inverse of the slope in 
log-linear scale) stay practically unchanged. Similarly, the 
vertical structure of the disk does not show significant evo- 
lution (top right panel) , except a moderate amount of flar- 
ing in the outer disk, which are due to spiral instabilities 
induced by swing-amplified Poisson noise in the disk. The 
scale-heights, measured at R = 2ARd, change from zo= 0.35 
kpc to 0.41 kpc in the "z=0" configuration and from 0.17 
kpc to 0.24 kpc in the "z=l" one. The disks are also slightly 
slowed down (middle left panel) , while the velocity disper- 
sions show an increase of ~ 5 km/s in the first Gyr in the 
fixed halo, and a total of < 10 km/s after 5 Gyr of evolution 
in the live potential (middle-right and bottom panels). The 
velocity ellipsoid of the disk (also measured at R = 2ARd) 
increases from (<7R,<70,cr,z)=(28,2O,17.5) km/s to (35,28,24) 
km/s in the "z=0" configuration and from (25,18,14) km/s 
to (32,28,22) km/s in the "z=l" one. 
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Figure A2. Evolution of the structural and kinematic proper- 
ties of the isolated host disk, for the "z=0" configuration. The 
solid lines show the initial conditions. The dotted lines show the 
disk after 1 Gyr of evolution within the fixed halo potential. The 
dashed lines show the disk's evolution after 5 Gyr within the N- 
body halo. The latter model is used as the control case in the 
comparison to the minor merger simulations. A similar behaviour 
is obtained for the disk in the "z=l" configuration. 



servations of energy and angular momentum that are better 
than 0.1% over 9 Gyr of evolution for our main disk galaxy 
configured at "z=0". 



A3 Evolution of Isolated Host 

To test the stability of the host galaxy, this is simulated in 
isolation, i.e. in the absence of any external perturbation. 

We first relax the disk component within a "rigid ver- 
sion" of the halo potential (which mimics the N-body halo 
described in ^Al.l|) for a few rotational periods (normally 
1 Gyr). Once the disk component is relaxed (i.e., there are 
no further changes in either its morphological or kinematical 
properties) the "rigid" halo is simply replaced by its A^-body 
("live") counterpart, and evolved for additional 5 Gyr in iso- 
lation for the configuration at "z=0" , and during 4 Gyr for 
that at "z=l". As described in H3.ll these time windows are 
enough to study the merger events that are of interest to us. 
Strong bar instabilities appear in the host galaxies evolved 



